#include "eb.h"

#ifdef ENABLE_EMPIRICAL_BAYES

#include "free.h"
#include "utilities.h"
#include "optimiz.h"
#include "pars.h"
#include "spr.h"
#include "lk.h"

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_rng.h>

//
// post-condition: writes an _eb_ file with several EB traces.
//
void Empirical_Bayes(t_tree* tree)
{
    int CHAIN_THINNING = 1;
    int DISPLAY_TIME_REMAINING_PERIOD = 10;
    int i;

    FILE* outfile = fopen(tree->io->out_eb_file, "w");

    if(!outfile)
    {
            fprintf(stderr, "ERROR: cant open output file for EB: %s\n", tree->io->out_eb_file);
            exit(1);
    }

    /* create random number generator */
    gsl_rng* rng = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(rng, (int)time(NULL));

    PhyML_Printf("\n. Starting empirical Bayes MCMC for %d generations.\n", tree->io->eb_ngens);
    PhyML_Printf("\n. Sampled trees during MCMC will be written to %s\n\n", tree->io->out_eb_file);

    /* Record clock time, to estimate time-to-completion */
    time_t start = time(NULL);
    time_t now;

    /* create storage for trees and models */
#ifndef ENABLE_MBL
	phydbl *best_bl = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
#else
	phydbl *best_bl = (phydbl *)mCalloc((2*tree->n_otu-3)*n_l,sizeof(phydbl));
#endif
	//PhyML_Printf("\n. debug eb.c 47 here");
	Record_Br_Len(best_bl,tree);
	//PhyML_Printf("\n. debug eb.c 49 here");
	t_tree *best_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
	//PhyML_Printf("\n. debug eb.c 51 here");
	Copy_Tree(tree, best_tree);
	//PhyML_Printf("\n. debug eb.c 53 here");
	model *best_model = Copy_Model(tree->mod);
	//PhyML_Printf("\n. debug eb.c 55 here");

	t_tree *last_tree = Make_Tree_From_Scratch(tree->n_otu,tree->data);
	model *last_model  = Copy_Model(tree->mod);
#ifndef ENABLE_MBL
	phydbl *last_bl = (phydbl *)mCalloc(2*tree->n_otu-3,sizeof(phydbl));
#else
	phydbl *last_bl = (phydbl *)mCalloc((2*tree->n_otu-3)*n_l,sizeof(phydbl));
#endif
	Record_Br_Len(last_bl,tree);
	Copy_Tree(tree, last_tree);

	//PhyML_Printf("\n. debug eb.c 58 here");

    double lnL_best = tree->c_lnL;
    double lnL_current = tree->c_lnL;


    /* MCMC routine */
    for(i=0; i<tree->io->eb_ngens; i++)
    {
            /* modify current tree */
            Random_Spr(tree);

            /* optimize current tree */
            Optimiz_All_Free_Param(tree,0);

            /* check acceptance */
            //printf("LR %f/%f [%f]\n", tree->c_lnL, lnL_current, exp(tree->c_lnL - lnL_current));

            if(tree->c_lnL >= lnL_current || gsl_rng_uniform(rng) < exp(tree->c_lnL - lnL_current))
            {
                    Copy_Tree(tree, last_tree);
                    Record_Model(tree->mod, last_model);
                    lnL_current = tree->c_lnL;
                    Record_Br_Len(last_bl, tree);
                    //Prepare_Tree_For_Lk(last_tree);


                    /* check if better than best so far */
                    if(tree->c_lnL > lnL_best)
                    {
                            Copy_Tree(tree, best_tree);
                            Record_Model(tree->mod, best_model);
                            lnL_best = tree->c_lnL;
                            Record_Br_Len(best_bl, tree);
                            //Prepare_Tree_For_Lk(best_tree);
                    }
            }
            else
            {
                    /* copy old tree back to current tree */
                    Copy_Tree(last_tree, tree);
                    Record_Model(last_model, tree->mod);
                    tree->c_lnL = lnL_current;
                    Restore_Br_Len(last_bl, tree);
                    //Prepare_Tree_For_Lk(tree);
            }

            /* print sampled tree */
            if(i % CHAIN_THINNING == 0)
            {
                    fprintf(outfile, "[%f] ", tree->c_lnL);
#ifndef ENABLE_MBL
                    Print_Tree(outfile, tree);
#else
                    Print_Tree(outfile, tree, -1);
#endif
            }

            /* print estimated remaining time*/
            if (i%DISPLAY_TIME_REMAINING_PERIOD == 0)
            {       now = time(NULL);
                    Print_time_remaining(now, start, i, tree->io->eb_ngens);
            }
    }

    /* copy best tree found back into tree */
    Copy_Tree(best_tree, tree);
    Record_Model(best_model, tree->mod);
    Prepare_Tree_For_Lk(tree);
    Restore_Br_Len(best_bl, tree);
    tree->c_lnL = lnL_best;

    /* clean up and done */
    fclose(outfile);
    gsl_rng_free(rng);
    PhyML_Printf("\n. Done with empirical Bayes MCMC.\n");

    Free_Model(best_model);
	Free_Tree(best_tree);
	Free_Model(last_model);
	Free_Tree(last_tree);
	Free(best_bl);
	Free(last_bl);
}

// Calculate posterior probabilities of clades,
// given the *.eb file written by the method named "Empricial_Bayes"
void PostProb_From_Tree(t_tree *tree, calign *alldata, model *mod, option *io)
{
	tree->both_sides = 1;
	Lk(tree);

	tree->print_alrt_val = 0;
	tree->print_boot_val = 0;
	tree->print_pp_val	 = 1;


	FILE* ebfile = fopen(tree->io->out_eb_file, "r");
	PhyML_Printf("\n. Calculating posterior probabilities of clades, using the MCMC samples in %s\n", tree->io->out_eb_file);

    if(!ebfile)
    {
            fprintf(stderr, "ERROR: cant open the file: %s\n",
            		tree->io->out_eb_file);
            exit(1);
    }

	int i,j,k;
	int count_lines = 0; // how many lines (i.e. samples) have we seen?
	int MAX_PARTS = 10 * (2 * (int)tree->n_otu - 2);
	struct __Node ***parts = (t_node ***)mCalloc(MAX_PARTS,sizeof(t_node **)); // parts[i] = a list_of_reachable_tips for partition i

	For(i, MAX_PARTS)
	{
		parts[i] = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *));
	}
	int *size_of_parts; // how big is each partition?
	size_of_parts = (int *)mCalloc(MAX_PARTS,sizeof(int));
	For(i,MAX_PARTS) size_of_parts[i] = 0;
	int *freq_of_parts; // how often have we seen each partition?
	freq_of_parts = (int *)mCalloc(MAX_PARTS,sizeof(int));
	For(i,MAX_PARTS) freq_of_parts[i] = 0;
	int count_parts = 0; // how many partitions exist?

	/*
	 * 1. Parse the MCMC samples and count the frequency of topological partitions among samples
	 */
	t_tree* sampled_tree = NULL;
	char line[T_MAX_LINE];
	while( fgets( line, T_MAX_LINE, ebfile) != NULL ){
		char *tokens = strtok( line, " ;");
		tokens = strtok( NULL, " ;");
		//sampled_tree = Restore_Tree_From_String(tokens, tree);
		sampled_tree = Read_Tree(tokens); // Read the Newick-formatted tree from the *.eb file and create an t_tree structure. . .
		sampled_tree->mod         = mod;
		sampled_tree->io          = io;
		sampled_tree->data        = alldata;
		sampled_tree->both_sides  = tree->both_sides;
		sampled_tree->n_pattern   = tree->n_pattern;
		Order_Tree_CSeq(tree,alldata);
		Fill_Dir_Table(sampled_tree);
		Update_Dirs(sampled_tree);

		count_lines++;

		if ( (count_lines != 0) && (count_lines%1000 == 0) )
        {       PhyML_Printf(" . [ MCMC sample %d ]\n", count_lines);
        }
		PhyML_Printf(".");

		For(j,2*sampled_tree->n_otu-3) // for every edge in the tree...
		{
			/*
			 * Does edge j connect to a terminal taxa?
			 * If so, then this partition will *always* have posterior probability = 1.0.
			 * Therefore, let's not bother computing the PP for this edge.
			 */
			if (sampled_tree->t_edges[j]->rght->tax || sampled_tree->t_edges[j]->left->tax)
			{	continue;
			}

			/*
			 * Here we find that largest partition attached to edge j.
			 * The tips in the largest partition will be saved in tiplist.
			 */
			t_node **tiplist;
			int proposed_part_size = 0;
			For(k,3)
			{
				if (sampled_tree->t_edges[j]->left->bip_size[k] > proposed_part_size)
				{	tiplist = sampled_tree->t_edges[j]->left->bip_node[k];
					proposed_part_size = sampled_tree->t_edges[j]->left->bip_size[k];
				}
				if (sampled_tree->t_edges[j]->rght->bip_size[k] > proposed_part_size)
				{	tiplist = sampled_tree->t_edges[j]->rght->bip_node[k];
					proposed_part_size = sampled_tree->t_edges[j]->rght->bip_size[k];
				}
			}

			// Now we compare our proposed partition to all other partitions found so far.
			// If our proposed partition is not unique, then we increment the counter for its match.
			// Else, we add the proposed partition to our list of partitions.
			int found_parti = 0;
			For(i,count_parts) // for each partition...
			{	// do partition i and subtree k have the same number of descendants?
				if (size_of_parts[i] == proposed_part_size)
				{	int score = Compare_List_Of_Reachable_Tips_version2(parts[i],
							size_of_parts[i],
							tiplist,
							proposed_part_size);

					// do partition i and subtree k contain the exact same set of descendant taxa?
					if(score == size_of_parts[i])
					{
						// ... then partition i is present within the sampled_tree.
						// Increment the observed frequency for partition i:
						found_parti = 1;
						freq_of_parts[i]++;
						break;
					}
					// a case we must deal with: if the partition size = 1/2 the number of terminal taxa, and if
					// the match score = 0, then we found the taxa in tiplist are the inverse of the taxa
					// in partition i.  In other words, we found the partition in sampled_tree.
					if (score == 0 && proposed_part_size == sampled_tree->n_otu * 0.5)
					{
						found_parti = 1;
						freq_of_parts[i]++;
						break;
					}
				}
			}
			if (found_parti == 0) // if the partition of node j has not been observed in previous samples...
			{
				// ... then add the partition to our collection of observed partitions.
				For(i, proposed_part_size)
				{
					parts[count_parts][i] = (t_node *)Make_Node_Light( tiplist[i]->num );
					strcpy(parts[count_parts][i]->name, tiplist[i]->name);
				}

				size_of_parts[count_parts] = proposed_part_size;
				freq_of_parts[count_parts]++; // ... and increment the observed frequency of this partition.
				count_parts++;
			}
		}
		if(sampled_tree->mat) Free_Mat(sampled_tree->mat);
		Free_Tree(sampled_tree);

	}
	PhyML_Printf("\n");
	fclose(ebfile);

	/*
	 * 2. Scale the raw frequency counts to posterior probabilities.
	 */
	double *pp_of_parts; // how often have we seen each partition?
	pp_of_parts = (double *)mCalloc(count_parts,sizeof(double));
	For(i,count_parts) pp_of_parts[i] = 0.0;
	For(i,count_parts)
	{	pp_of_parts[i] = (double)freq_of_parts[i] / (double)count_lines;
	}

	/*
	 * 3. Write posterior probabilities onto the branches of the given tree
	 */
	For(j,2*tree->n_otu-3) // for every e dge in the tree...
	{
		/*
		 * Does edge j connect to a terminal taxa?
		 * If so, then this partition will *always* have posterior probability = 1.0.
		 * Therefore, let's not bother computing the PP for this edge.
		 */
		if (tree->t_edges[j]->rght->tax || tree->t_edges[j]->left->tax)
		{	tree->t_edges[j]->post_prob = 1.0;
		}

		/*
		 * Determine which partition matches edge j
		 */
		t_node **tiplist;
		int proposed_part_size = 0;
		For(k,3)
		{
			if (tree->t_edges[j]->left->bip_size[k] > proposed_part_size)
			{	tiplist = tree->t_edges[j]->left->bip_node[k];
				proposed_part_size = tree->t_edges[j]->left->bip_size[k];
			}
			if (tree->t_edges[j]->rght->bip_size[k] > proposed_part_size)
			{	tiplist = tree->t_edges[j]->rght->bip_node[k];
				proposed_part_size = tree->t_edges[j]->rght->bip_size[k];
			}
		}
		For(i,count_parts)
		{	if (size_of_parts[i] == proposed_part_size)
			{
				int score = Compare_List_Of_Reachable_Tips_version2(parts[i],
						size_of_parts[i],
						tiplist,
						proposed_part_size);

				if( (score == size_of_parts[i]) ||
						(score == 0 && proposed_part_size == tree->n_otu * 0.5) )
				{	tree->t_edges[j]->post_prob = pp_of_parts[i];
					break;
				}
			}
		}
	}


	/*
	 * 4. Write a partition table
	 *
	 * The relationship between taxa names and numbers will accord with the order of tree->noeud[i]
	 */
	FILE* partsfile = fopen(tree->io->out_parts_file, "w");
    if(!partsfile)
    {
            fprintf(stderr, "ERROR: cant open the file: %s\n",
            		tree->io->out_parts_file);
            exit(1);
    }


	// Print a list of all taxa numbers
	For(i, tree->n_otu)
	{	fprintf(partsfile, "Taxa #%d = %s\n", i, tree->noeud[i]->name);
	}

	fprintf(partsfile, "\nPartitions:\n");
	fprintf(partsfile, "\t* = in partition\n\t- = excluded from partition\n");
	fprintf(partsfile, "\tFreq. = the number of samples in which this partition appears.\n");
	fprintf(partsfile, "\tPP = the posterior probability of this partition.\n");
	fprintf(partsfile, "\n");

	For(i,count_parts)
	{
		PhyML_Printf(" . debug, partition %d\n", i);

		fprintf(partsfile, "[p%d]\t", i);
		int *part = (int *)mCalloc(tree->n_otu,sizeof(int));
		For(j,tree->n_otu)
		{	part[j] = 0;
		}
		For(j,size_of_parts[i])
		{	//PhyML_Printf("%s|%d ",   parts[i][j]->name, parts[i][j]->num);

			// which 'num' is this node?
			int num = -1;
			For(k, tree->n_otu)
			{
				if (0 == strcmp(parts[i][j]->name, tree->noeud[k]->name) )
				{
					num = k;
					break;
				}
			}
			part[ num ] = 1; // mark node #num as being in this partition
		}

		For(j,tree->n_otu)
		{
			if (part[j] == 0)
			{	fprintf(partsfile, "*");
			}
			else if (part[j] == 1 )
			{	fprintf(partsfile, "-");
			}
		}
		PhyML_Printf("\t[Freq. = %d]\t [PP = %f]\n", freq_of_parts[i], pp_of_parts[i]);
		fprintf(partsfile, "\t[Freq. = %d]\t [PP = %f]\n", freq_of_parts[i], pp_of_parts[i]);
		Free(part);

	}

	fprintf(partsfile, "\n");
	fclose(partsfile);
	For(i, MAX_PARTS)
	{
		For(j, size_of_parts[i])
		{
			Free_Node(parts[i][j]);
		}
	}
	Free(parts);

	Free(size_of_parts);
	Free(freq_of_parts);
	Free(pp_of_parts);
}
/*
	Free(s_tree);
	s_tree = Write_Tree(tree,1);

	Free_Spr_List(tree);
	Free_One_Spr(tree->best_spr);
	Free_Triplet(tree->triplet_struct);
	Free_Tree_Pars(tree);
	Free_Tree_Lk(tree);
	Free_Tree(tree);

	return s_tree;
}
*/

// Calculate posterior probabilities of clades,
// given the *.eb file written by the method named "Empricial_Bayes"
char *PostProb_From_String(char *s_tree, calign *alldata, model *mod, option *io)
{
	t_tree *tree;
	tree = Read_Tree(s_tree);

	if(!tree)
	{
		PhyML_Printf("\n. Err in file %s at line %d\n",__FILE__,__LINE__);
		Warn_And_Exit("");
	}

	tree->mod         = mod;
	tree->io          = io;
	tree->data        = alldata;
	tree->both_sides  = 1;
	tree->n_pattern   = tree->data->crunch_len;

	Order_Tree_CSeq(tree,alldata);
	if(tree->mod->s_opt->random_input_tree) Random_Tree(tree);
	Fill_Dir_Table(tree);
	Update_Dirs(tree);
	Make_Tree_4_Pars(tree,alldata,alldata->init_len);
	Make_Tree_4_Lk(tree,alldata,alldata->init_len);
	tree->triplet_struct = Make_Triplet_Struct(mod);
	Br_Len_Not_Involving_Invar(tree);
	Make_Spr_List(tree);
	Make_Best_Spr(tree);

	tree->both_sides = 1;
	Lk(tree);

	tree->print_alrt_val = 0;
	tree->print_boot_val = 0;
	tree->print_pp_val	 = 1;

	FILE* ebfile = NULL;

	ebfile = fopen(tree->io->out_eb_file, "r");
	PhyML_Printf("\n. Calculating posterior probabilities of clades, using the MCMC samples in %s\n", tree->io->out_eb_file);

	if(!ebfile)
	{
        fprintf(stderr, "ERROR: cant open the file: %s\n",
        		tree->io->out_eb_file);
		exit(1);
	}

	int i,j,k;
	int count_lines = 0; // how many lines (i.e. samples) have we seen?
	int MAX_PARTS = 10 * (2 * (int)tree->n_otu - 2);
	struct __Node ***parts = (t_node ***)mCalloc(MAX_PARTS,sizeof(t_node **)); // parts[i] = a list_of_reachable_tips for partition i

	For(i, MAX_PARTS)
	{
		parts[i] = (t_node **)mCalloc(tree->n_otu,sizeof(t_node *));
	}
	int *size_of_parts; // how big is each partition?
	size_of_parts = (int *)mCalloc(MAX_PARTS,sizeof(int));
	For(i,MAX_PARTS) size_of_parts[i] = 0;
	int *freq_of_parts; // how often have we seen each partition?
	freq_of_parts = (int *)mCalloc(MAX_PARTS,sizeof(int));
	For(i,MAX_PARTS) freq_of_parts[i] = 0;
	int count_parts = 0; // how many partitions exist?

	/*
	 * 1. Parse the MCMC samples and count the frequency of topological partitions among samples
	 */
	t_tree* sampled_tree = NULL;
	char line[T_MAX_LINE];
	while( fgets( line, T_MAX_LINE, ebfile) != NULL ){
		char *tokens = strtok( line, " ;");
		//PhyML_Printf("\n. eb.c 581, line=%s", tokens);
		tokens = strtok( NULL, " ;");
		sampled_tree = Read_Tree(tokens); // Read the Newick-formatted tree from the *.eb file and create an t_tree structure. . .
		sampled_tree->mod         = tree->mod;
		sampled_tree->io          = tree->io;
		sampled_tree->data        = tree->data;
		sampled_tree->both_sides  = tree->both_sides;
		sampled_tree->n_pattern   = tree->n_pattern;
		Fill_Dir_Table(sampled_tree);
		Update_Dirs(sampled_tree);

		count_lines++;

		if ( (count_lines != 0) && (count_lines%1000 == 0) )
        {       PhyML_Printf(" . [ MCMC sample %d ]\n", count_lines);
        }
		//PhyML_Printf("*");

		For(j,2*sampled_tree->n_otu-3) // for every edge in the tree...
		{
			/*
			 * Does edge j connect to a terminal taxa?
			 * If so, then this partition will *always* have posterior probability = 1.0.
			 * Therefore, let's not bother computing the PP for this edge.
			 */
			if (sampled_tree->t_edges[j]->rght->tax || sampled_tree->t_edges[j]->left->tax)
			{	continue;
			}

			//PhyML_Printf("\n. length of sampled edge = %f", sampled_tree->t_edges[j]->l);

			/*
			 * Here we find that largest partition attached to edge j.
			 * The tips in the largest partition will be saved in tiplist.
			 */
			t_node **tiplist;
			int proposed_part_size = sampled_tree->n_otu;
			For(k,3)
			{
				//PhyML_Printf("\n. eb.c 563 -, edge %d, k=%d, left bip size = %d, right bip size = %d", j, k, sampled_tree->t_edges[j]->left->bip_size[k], sampled_tree->t_edges[j]->rght->bip_size[k]);

				// This line is a wierd hack: VHS
				// Unresolved branches appear to evade my previous test (see: Does edge j connect to a terminal taxa?)
				// Here we ensure the partition is not a terminal edge
				if (sampled_tree->t_edges[j]->rght->bip_size[k] == sampled_tree->n_otu-1 || sampled_tree->t_edges[j]->left->bip_size[k] == sampled_tree->n_otu-1)
				{	continue;
				}

				if (sampled_tree->t_edges[j]->left->bip_size[k] < proposed_part_size)
				{	tiplist = sampled_tree->t_edges[j]->left->bip_node[k];
					proposed_part_size = sampled_tree->t_edges[j]->left->bip_size[k];
				}
				if (sampled_tree->t_edges[j]->rght->bip_size[k] < proposed_part_size)
				{	tiplist = sampled_tree->t_edges[j]->rght->bip_node[k];
					proposed_part_size = sampled_tree->t_edges[j]->rght->bip_size[k];
				}
			}

			//PhyML_Printf("\n. proposed_part_size = %d", proposed_part_size);

			// for debug:
			/*
			int qqq;
			PhyML_Printf("\n. tiplist:");
			for (qqq=0; qqq< proposed_part_size; qqq++)
			{
				PhyML_Printf("\n. %s", tiplist[qqq]->name);
			}
			PhyML_Printf("\n. eb.c 565 - proposed_part_size = %d", proposed_part_size);
			PhyML_Printf("\n. eb.c 566 - count_parts = %d", count_parts);
			 */

			// Now we compare our proposed partition to all other partitions found so far.
			// If our proposed partition is not unique, then we increment the counter for its match.
			int found_parti = 0;
			For(i,count_parts) // for each partition...
			{	// do partition i and subtree k have the same number of descendants?
				if (size_of_parts[i] == proposed_part_size)
				{	int score = Compare_List_Of_Reachable_Tips_version2(parts[i],
							size_of_parts[i],
							tiplist,
							proposed_part_size);

					// do partition i and subtree k contain the exact same set of descendant taxa?
					if(score == size_of_parts[i])
					{
						// ... then partition i is present within the sampled_tree.
						// Increment the observed frequency for partition i:
						found_parti = 1;
						freq_of_parts[i]++;
						break;
					}
					// a case we must deal with: if the partition size = 1/2 the number of terminal taxa, and if
					// the match score == 0, then we found the taxa in tiplist are the inverse of the taxa
					// in partition i.  In other words, we found the partition in sampled_tree.
					else if (score == 0 && proposed_part_size == sampled_tree->n_otu * 0.5)
					{
						found_parti = 1;
						freq_of_parts[i]++;
						break;
					}
				}
			}
			// If we didn't find the partition in the above For block,
			// then here we add the proposed partition to our list of partitions.
			//PhyML_Printf("\n. eb.c 598");
			if (found_parti == 0) // if the partition of node j has not been observed in previous samples...
			{
				// ... then add the partition to our collection of observed partitions.
				For(i, proposed_part_size)
				{
					parts[count_parts][i] = (t_node *)Make_Node_Light( tiplist[i]->num );
					parts[count_parts][i]->name = (char *)mCalloc((int)strlen(tiplist[i]->name)+1,sizeof(char));
					strcpy(parts[count_parts][i]->name, tiplist[i]->name);
				}

				size_of_parts[count_parts] = proposed_part_size;
				freq_of_parts[count_parts]++; // ... and increment the observed frequency of this partition.
				count_parts++;
			}
			//PhyML_Printf("\n. eb.c 612");
		}
		if(tree->mat) Free_Mat(sampled_tree->mat);
		Free_Tree(sampled_tree);

	}
	PhyML_Printf("\n");
	fclose(ebfile);

	/*
	 * 2. Scale the raw frequency counts to posterior probabilities.
	 */
	double *pp_of_parts; // how often have we seen each partition?
	pp_of_parts = (double *)mCalloc(count_parts,sizeof(double));
	For(i,count_parts) pp_of_parts[i] = 0.0;
	For(i,count_parts)
	{	pp_of_parts[i] = (double)freq_of_parts[i] / (double)count_lines;
	}

	/*
	 * 3. Write posterior probabilities onto the branches of the given tree
	 */
	For(j,2*tree->n_otu-3) // for every e dge in the tree...
	{
		/*
		 * Does edge j connect to a terminal taxa?
		 * If so, then this partition will *always* have posterior probability = 1.0.
		 * Therefore, let's not bother computing the PP for this edge.
		 */
		if (tree->t_edges[j]->rght->tax || tree->t_edges[j]->left->tax)
		{	tree->t_edges[j]->post_prob = 1.0;
		}

		/*
		 * Determine which partition matches edge j
		 */
		t_node **tiplist;
		int proposed_part_size = tree->n_otu;
		For(k,3)
		{
			if (tree->t_edges[j]->rght->bip_size[k] == tree->n_otu-1 || tree->t_edges[j]->left->bip_size[k] == tree->n_otu-1)
			{	continue;
			}

			if (tree->t_edges[j]->left->bip_size[k] < proposed_part_size)
			{	tiplist = tree->t_edges[j]->left->bip_node[k];
				proposed_part_size = tree->t_edges[j]->left->bip_size[k];
			}
			if (tree->t_edges[j]->rght->bip_size[k] < proposed_part_size)
			{	tiplist = tree->t_edges[j]->rght->bip_node[k];
				proposed_part_size = tree->t_edges[j]->rght->bip_size[k];
			}
		}
		For(i,count_parts)
		{	if (size_of_parts[i] == proposed_part_size)
			{
				int score = Compare_List_Of_Reachable_Tips_version2(parts[i],
						size_of_parts[i],
						tiplist,
						proposed_part_size);

				if( (score == size_of_parts[i]) ||
						(score == 0 && proposed_part_size == tree->n_otu * 0.5) )
				{	tree->t_edges[j]->post_prob = pp_of_parts[i];
					break;
				}
			}
		}
	}


	/*
	 * 4. Write a partition table
	 *
	 * The relationship between taxa names and numbers will accord with the order of tree->noeud[i]
	 */
    PhyML_Printf("\n. Printing MCMC sampled tree partitions to '%s'...\n\n", tree->io->out_parts_file);
	FILE* partsfile = fopen(tree->io->out_parts_file, "w");
    if(!partsfile)
    {
            fprintf(stderr, "ERROR: cant open the file: %s\n",
            		tree->io->out_parts_file);
            exit(1);
    }

	// Print a list of all taxa numbers
	For(i, tree->n_otu)
	{	fprintf(partsfile, "Taxa #%d = %s\n", i, tree->noeud[i]->name);
	}

	fprintf(partsfile, "\nPartitions:\n");
	fprintf(partsfile, "\t* = in partition\n\t- = excluded from partition\n");
	fprintf(partsfile, "\tFreq. = the number of samples in which this partition appears.\n");
	fprintf(partsfile, "\tPP = the posterior probability of this partition.\n");
	fprintf(partsfile, "\n");

	For(i,count_parts)
	{
		//PhyML_Printf(" . debug, partition %d\n", i);

		fprintf(partsfile, "[p%d]\t", i);
		int *part = (int *)mCalloc(tree->n_otu,sizeof(int));
		For(j,tree->n_otu)
		{	part[j] = 0;
		}
		For(j,size_of_parts[i])
		{	//PhyML_Printf("%s|%d ",   parts[i][j]->name, parts[i][j]->num);

			// which 'num' is this node?
			int num = -1;
			For(k, tree->n_otu)
			{
				if (0 == strcmp(parts[i][j]->name, tree->noeud[k]->name) )
				{
					num = k;
					break;
				}
			}
			part[ num ] = 1; // mark node #num as being in this partition
		}

		For(j,tree->n_otu)
		{
			if (part[j] == 0)
			{	fprintf(partsfile, "*");
			}
			else if (part[j] == 1 )
			{	fprintf(partsfile, "-");
			}
		}
		fprintf(partsfile, "\t[Freq. = %d]\t [PP = %f]\n", freq_of_parts[i], pp_of_parts[i]);
		Free(part);

	}

	fprintf(partsfile, "\n");
	fclose(partsfile);
	For(i, MAX_PARTS)
	{
		For(j, size_of_parts[i])
		{
			Free_Node(parts[i][j]);
		}
	}
	Free(parts);

	Free(size_of_parts);
	Free(freq_of_parts);
	Free(pp_of_parts);

	Free(s_tree);
#ifndef ENABLE_MBL
	s_tree = Write_Tree(tree);
#else
	s_tree = Write_Tree(tree,-1);
#endif

	Free_Spr_List(tree);
	Free_One_Spr(tree->best_spr);
	//Free_Triplet(tree->triplet_struct);
	Free_Tree_Pars(tree);
	Free_Tree_Lk(tree);
	Free_Tree(tree);

	return s_tree;
}

/*********************************************************/
// now = the time_t object recorded now
// start = the time_t object recorded when we started
// i = the iteration number of the current generation, where 0 <= i < total
// total = the total number of generations we want to complete
void Print_time_remaining(time_t now, time_t start, int i, int total)
{
	// etc = "estimated time to completion" in seconds
	long delta = (long)difftime(now,start);
	phydbl secs_per_rep = delta / (phydbl)(i+1);
	int remaining_reps = total - (i + 1);
	int seconds = remaining_reps * secs_per_rep;

	int SECONDS_IN_DAY = 86400;
	int SECONDS_IN_HOUR = 3600;
	int SECONDS_IN_MINUTE = 60;

	int days = seconds / SECONDS_IN_DAY;
	seconds = seconds - (days * SECONDS_IN_DAY);
	int hours = seconds / SECONDS_IN_HOUR;
	seconds = seconds - (hours * SECONDS_IN_HOUR);
	int minutes = seconds / SECONDS_IN_MINUTE;
	seconds = seconds - (minutes * SECONDS_IN_MINUTE);

    if (days > 0)
    {
    	PhyML_Printf("\n . EB Est. time remaining = %dd %dh %dm %lds\n", days, hours, minutes, seconds );
    }
    else
    {
    	PhyML_Printf("\n . EB Est. time remaining = %dh %dm %lds \n", hours, minutes, seconds );
    }
}

/*
 * Victor: I  wrote this version as an alternative to Compare_List_Of_Reachable_Tips.
 * When I use Compare_List_Of_Reachable_Tips in our empirical Bayesian methods
 * (see Calculate_PP), Compare_List_Of_Reachable_Tips returns incorrect match counts.
 * Perhaps my lists are formatted incorrectly?
 * In this alternative version, instead of comparing tip objects, I compare tip names.
 */
int Compare_List_Of_Reachable_Tips_version2(t_node **list1, int size_list1, t_node **list2, int size_list2)
{
	//PhyML_Printf(". debug utilities.c 7442 entered Compare_list\n");
	int i,j,n_matches;

	n_matches = 0;
	For(i,size_list1)
	{
		For(j,size_list2)
		{
			/* The following line is different from the original version of Compare_List_Of_Reachable_Tips.
			 * By the way, strcmp doesn't return typical 0 = false, 1 = true values.
			 * Instead, 0 = the strings are identical.  For more information, see this reference:
			 * http://www.cplusplus.com/reference/clibrary/cstring/strcmp/
			 */

			if (0 == strcmp(list1[i]->name, list2[j]->name) )
			{
				n_matches++;
			}

			//PhyML_Printf(" %d = strcmp( %s, %s )\n", strcmp(list1[i]->name, list2[j]->name),  list1[i]->name, list2[j]->name);
		}
	}
	//PhyML_Printf(". debug utilities.c 7646 leaving Compare_list\n");
	return n_matches;
}

/*********************************************************/

void Random_Spr(t_tree *tree)
{
	// pick a random branch:
	unsigned int iseed = (unsigned int)time(NULL);
	srand(iseed);
	int randb = rand() % (2*tree->n_otu-3); // this is the index of a random branch

	t_edge *b = tree->t_edges[randb];
	t_node *n = b->left;
	if (b->left->tax) n = b->rght;
	// VHS: to-do:
	// NOTE: We probably simplify this next step, by not testing all possible targets,
	// but rather picking just one random target.
	Test_All_Spr_Targets(b,n,tree);

	// Make a random SPR target from the list of all possible targets that were
	// computed by Test_All_Spr_Targets.
	int found_good_move = 0;
	int rand_movei = 0;
	for (rand_movei = 0; rand_movei < tree->size_spr_list; rand_movei++)
	{

		//char* s_tree = (char *)Write_Tree(tree);
		//PhyML_Printf("\n. %s", s_tree);

		//int rand_movei = rand() % tree->size_spr_list;
		//int rand_movei = 0;
		spr *move = tree->spr_list[rand_movei];
		if (move->n_link != move->n_opp_to_link)
		{
			found_good_move = 1;
			//PhyML_Printf("\n. eb.c using random move %d, n_link=%d, n_opp_to_link=%d", rand_movei, move->n_link, move->n_opp_to_link);
			t_edge *target,*residual;
			Prune_Subtree(move->n_link,move->n_opp_to_link,&target,&residual,tree);
			//PhyML_Printf("\n. eb.c 874");
			Graft_Subtree(move->b_target,move->n_link,residual,tree);
			//PhyML_Printf("\n. eb.c 876");
			tree->both_sides = 1;
			//PhyML_Printf("\n. eb.c 878");
			Pars(tree);
			tree->best_pars = tree->c_pars;
		}
		if (found_good_move == 1)
		{
			break;
		}
	}
    Reset_Spr_List(tree);
}


#endif



